iSSF Vegetation Model Plotting

Author

Scott Forrest and Kyana Pike

Published

May 20, 2025

Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain vegetation categories as a covariate, with interactions for movement parameters.

Load packages

Code
library(tidyverse)
library(TwoStepCLogit) # to make population estimates
library(sjPlot) # to plot model estimates
library(glmmTMB)# for the mixed models 
library(tictoc)
library(beepr)
library(lemon)
library(circular)
library(amt)

Import data

Code
ssf_dat <- read_csv("pheno_start_end_dat_ssf10rs_2024-03-14.csv")
New names:
Rows: 1062908 Columns: 42
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(4): pheno_start, pheno_end, season, day_period dbl (35): ...1, id, burst_,
x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr... lgl (1): case_ dttm (2):
t1_, t2_
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
Code
# ssf_dat <- read_csv("pheno_forest_start_end_dat_ssf10rs_2025-03-02.csv")

Prepare the data

Code
# scale the covariates
ssf_dat <- ssf_dat %>% mutate(srtm_start_scaled = scale(srtm_start),
                            srtm_end_scaled = scale(srtm_end),
                            water_dist_start_scaled = scale(water_dist_start),
                            water_dist_end_scaled = scale(water_dist_end),
                            ndvi_start_scaled = scale(ndvi_start),
                            ndvi_end_scaled = scale(ndvi_end),
                            nbr_start_scaled = scale(nbr_start),
                            nbr_end_scaled = scale(nbr_end),
                            sl_scaled = scale(sl_),
                            log_sl_scaled = scale(log_sl_),
                            cos_ta_scaled = scale(cos_ta_)
)

# get scaling factors
sl_scale_sd <- attr(ssf_dat$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_dat$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_dat$cos_ta_scaled, "scaled:scale")

# rename the pheno categories to set a different reference level
ssf_dat$pheno_end <- relevel(as.factor(ssf_dat$pheno_end), ref = "shrub_grass")
ssf_dat$pheno_start <- relevel(as.factor(ssf_dat$pheno_start), ref = "shrub_grass")

unique(ssf_dat$pheno_start)
[1] shrub_grass   dry_grassland tree_grass    bare_ground   open_forest  
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass
Code
unique(ssf_dat$pheno_end)
[1] shrub_grass   dry_grassland tree_grass    bare_ground   open_forest  
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass
Code
ssf_dat <- ssf_dat %>% filter(!pheno_start == "saltpan_mudlfat")

# separate the buffalo by season 
ssfdat_wet <- ssf_dat %>% filter(season=="wet") # wet season
ssfdat_dry <- ssf_dat %>% filter(season=="dry") # dry season only 

Plot the used-available distributions

Code
# checking for all individuals
buffalo_ids <- unique(ssf_dat$id)

# inspect the data for dry season 
ssfdat_dry %>% 
  group_by(case_, pheno_end) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season", subtitle = "All buffalo")+
  scale_fill_brewer(palette = "Paired", name="Case", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code
ggsave("outputs/plots/used-available_all_dry.png", device = "png", 
       width = 150, height = 90, units = "mm", scale = 1.5, dpi = 1000)

Inspect the data for each buffalo

Code
for(i in 1:length(buffalo_ids)) {
  print(ssfdat_dry %>% filter(id == buffalo_ids[i]) %>%
          group_by(case_, pheno_end) %>% 
          summarize(n = n()) %>% 
          mutate(prop = n / sum(n), 
                 label = paste0(round(prop * 100, 1), "%")) %>% 
          ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
          geom_col(position = position_dodge2()) +
          geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
          labs(x = "Land use class", y = "Proportion", fill = "case_")+
          ggtitle("Dry season", subtitle = paste0("Buffalo ", buffalo_ids[i])) +
          scale_fill_brewer(palette = "Paired", name="Case", 
                            breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
          theme_light() +
          theme(axis.text.x = element_text(angle = 20, hjust = 1))
  )
}
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code
ssfdat_wet %>% 
  group_by(case_, pheno_end) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "case_")+
  ggtitle("Wet season", subtitle = "All buffalo")+
  scale_fill_brewer(palette = "Paired", name="Case", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code
ggsave("outputs/plots/used-available_all_wet.png", device = "png", 
       width = 150, height = 90, units = "mm", scale = 1.5, dpi = 1000)

For each individual

Code
for(i in 1:length(buffalo_ids)) {
  print(ssfdat_wet %>% filter(id == buffalo_ids[i]) %>%
          group_by(case_, pheno_end) %>% 
          summarize(n = n()) %>% 
          mutate(prop = n / sum(n), 
                 label = paste0(round(prop * 100, 1), "%")) %>% 
          ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
          geom_col(position = position_dodge2()) +
          geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
          labs(x = "Land use class", y = "Proportion", fill = "case_")+
          ggtitle("Wet season", subtitle = paste0("Buffalo ", buffalo_ids[i]))+
          scale_fill_brewer(palette = "Paired", name="Case", 
                            breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
          theme_light() +
          theme(axis.text.x = element_text(angle = 20, hjust = 1))
  )
}
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Read in the saved model objects

Code
# buffalo.pheno.dry.mvmt <- readRDS("buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_dry_2025-05-20.rds")
# buffalo.pheno.wet.mvmt <- readRDS("buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_wet_2025-05-20.rds")

Model summaries

Code
# Summaries are also printed below
# summary(buffalo.pheno.dry.mvmt)
# summary(buffalo.pheno.wet.mvmt)

Model fitting - interacting with movement

Code
tic()
buffalo.tmp <- glmmTMB(case_ ~ 
                         pheno_end + 
                         sl_scaled +
                         log_sl_scaled +
                         cos_ta_scaled +
                         pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
                         (1 | step_id) + 
                         # (0 + pheno_end | id),
                         diag(0 + pheno_end | id),
                       family = poisson(), 
                       data = ssfdat_dry, 
                       doFit=FALSE,
                       control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))

#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA,1:(nvarparm-1))))

# to render the script we comment out the model fitting, as we are loading a saved model
buffalo.pheno.dry.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)

summary(buffalo.pheno.dry.mvmt) 
 Family: poisson  ( log )
Formula:          
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +  
    (1 | step_id) + diag(0 + pheno_end | id)
Data: ssfdat_dry

      AIC       BIC    logLik -2*log(L)  df.resid 
1160974.1 1161257.6 -580462.1 1160924.1    620298 

Random effects:

Conditional model:
 Groups  Name                   Variance  Std.Dev.  Corr                
 step_id (Intercept)            1.000e+06 1000.0000                     
 id      pheno_endshrub_grass   4.206e-02    0.2051                     
         pheno_endbare_ground   1.289e-02    0.1136 0.00                
         pheno_enddry_grassland 2.809e-02    0.1676 0.00 0.00           
         pheno_endopen_forest   1.762e-01    0.4197 0.00 0.00 0.00      
         pheno_endtree_grass    2.464e-02    0.1570 0.00 0.00 0.00 0.00 
Number of obs: 620323, groups:  step_id, 56393; id, 15

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                            -2.905077   4.211376  -0.690 0.490310
pheno_endbare_ground                    0.289401   0.064373   4.496 6.93e-06
pheno_enddry_grassland                  0.026845   0.071515   0.375 0.707379
pheno_endopen_forest                   -0.023060   0.123586  -0.187 0.851983
pheno_endtree_grass                    -0.011854   0.070324  -0.169 0.866140
sl_scaled                               0.017248   0.007990   2.159 0.030861
log_sl_scaled                          -0.003373   0.008243  -0.409 0.682367
cos_ta_scaled                          -0.001107   0.006586  -0.168 0.866556
sl_scaled:pheno_startbare_ground       -0.024199   0.017730  -1.365 0.172290
sl_scaled:pheno_startdry_grassland      0.058082   0.016770   3.464 0.000533
sl_scaled:pheno_startopen_forest       -0.033550   0.016937  -1.981 0.047607
sl_scaled:pheno_starttree_grass         0.018553   0.018620   0.996 0.319039
log_sl_scaled:pheno_startbare_ground   -0.040732   0.015809  -2.577 0.009981
log_sl_scaled:pheno_startdry_grassland -0.094216   0.017167  -5.488 4.06e-08
log_sl_scaled:pheno_startopen_forest    0.074716   0.017826   4.191 2.77e-05
log_sl_scaled:pheno_starttree_grass     0.013739   0.019710   0.697 0.485745
cos_ta_scaled:pheno_startbare_ground   -0.046518   0.012831  -3.625 0.000289
cos_ta_scaled:pheno_startdry_grassland  0.009236   0.013844   0.667 0.504667
cos_ta_scaled:pheno_startopen_forest   -0.038082   0.013482  -2.825 0.004732
cos_ta_scaled:pheno_starttree_grass    -0.035782   0.014887  -2.404 0.016233
                                          
(Intercept)                               
pheno_endbare_ground                   ***
pheno_enddry_grassland                    
pheno_endopen_forest                      
pheno_endtree_grass                       
sl_scaled                              *  
log_sl_scaled                             
cos_ta_scaled                             
sl_scaled:pheno_startbare_ground          
sl_scaled:pheno_startdry_grassland     ***
sl_scaled:pheno_startopen_forest       *  
sl_scaled:pheno_starttree_grass           
log_sl_scaled:pheno_startbare_ground   ** 
log_sl_scaled:pheno_startdry_grassland ***
log_sl_scaled:pheno_startopen_forest   ***
log_sl_scaled:pheno_starttree_grass       
cos_ta_scaled:pheno_startbare_ground   ***
cos_ta_scaled:pheno_startdry_grassland    
cos_ta_scaled:pheno_startopen_forest   ** 
cos_ta_scaled:pheno_starttree_grass    *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
toc()
130.739 sec elapsed
Code
beep(sound = 2)

# plot result as incidence rate ratios
# plot_model(buffalo.pheno.dry.mvmt)
plot_model(buffalo.pheno.dry.mvmt, transform = NULL)

Code
saveRDS(buffalo.pheno.dry.mvmt, file = paste0("buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_dry_", Sys.Date(), ".rds"))
paste0("Model location: buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_dry_", Sys.Date(), ".rds")
[1] "Model location: buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_dry_2025-05-20.rds"
Code
tic()
buffalo.tmp <- glmmTMB(case_ ~ 
                         pheno_end + 
                         sl_scaled +
                         log_sl_scaled +
                         cos_ta_scaled +
                         pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
                         (1 | step_id) + 
                         # (0 + pheno_end | id),
                         diag(0 + pheno_end | id),
                       family=poisson(), 
                       data = ssfdat_wet, 
                       doFit=FALSE,
                       control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))

#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA,1:(nvarparm-1))))

# to render the script we comment out the model fitting, as we are loading a saved model
buffalo.pheno.wet.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)

summary(buffalo.pheno.wet.mvmt)
 Family: poisson  ( log )
Formula:          
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +  
    (1 | step_id) + diag(0 + pheno_end | id)
Data: ssfdat_wet

      AIC       BIC    logLik -2*log(L)  df.resid 
 827189.3  827464.3 -413569.6  827139.3    442560 

Random effects:

Conditional model:
 Groups  Name                   Variance  Std.Dev.  Corr                
 step_id (Intercept)            1.000e+06 1.000e+03                     
 id      pheno_endshrub_grass   2.569e-02 1.603e-01                     
         pheno_endbare_ground   7.612e-02 2.759e-01 0.00                
         pheno_enddry_grassland 2.976e-05 5.455e-03 0.00 0.00           
         pheno_endopen_forest   7.819e-02 2.796e-01 0.00 0.00 0.00      
         pheno_endtree_grass    1.721e-02 1.312e-01 0.00 0.00 0.00 0.00 
Number of obs: 442585, groups:  step_id, 40235; id, 15

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                            -3.142132   4.985576  -0.630 0.528535
pheno_endbare_ground                    0.500199   0.087277   5.731 9.97e-09
pheno_enddry_grassland                  0.070347   0.048376   1.454 0.145903
pheno_endopen_forest                    0.019023   0.090207   0.211 0.832977
pheno_endtree_grass                     0.048338   0.060670   0.797 0.425597
sl_scaled                              -0.001803   0.008832  -0.204 0.838234
log_sl_scaled                           0.154034   0.011179  13.779  < 2e-16
cos_ta_scaled                           0.084257   0.007907  10.656  < 2e-16
sl_scaled:pheno_startbare_ground       -0.009350   0.018311  -0.511 0.609606
sl_scaled:pheno_startdry_grassland      0.007096   0.021411   0.331 0.740318
sl_scaled:pheno_startopen_forest       -0.090468   0.022801  -3.968 7.25e-05
sl_scaled:pheno_starttree_grass        -0.040275   0.025024  -1.609 0.107510
log_sl_scaled:pheno_startbare_ground   -0.074320   0.019119  -3.887 0.000101
log_sl_scaled:pheno_startdry_grassland -0.051027   0.025076  -2.035 0.041866
log_sl_scaled:pheno_startopen_forest   -0.114017   0.023049  -4.947 7.55e-07
log_sl_scaled:pheno_starttree_grass    -0.156630   0.024821  -6.310 2.78e-10
cos_ta_scaled:pheno_startbare_ground   -0.104274   0.014004  -7.446 9.63e-14
cos_ta_scaled:pheno_startdry_grassland -0.013069   0.017470  -0.748 0.454429
cos_ta_scaled:pheno_startopen_forest   -0.156955   0.017241  -9.104  < 2e-16
cos_ta_scaled:pheno_starttree_grass    -0.141402   0.018059  -7.830 4.89e-15
                                          
(Intercept)                               
pheno_endbare_ground                   ***
pheno_enddry_grassland                    
pheno_endopen_forest                      
pheno_endtree_grass                       
sl_scaled                                 
log_sl_scaled                          ***
cos_ta_scaled                          ***
sl_scaled:pheno_startbare_ground          
sl_scaled:pheno_startdry_grassland        
sl_scaled:pheno_startopen_forest       ***
sl_scaled:pheno_starttree_grass           
log_sl_scaled:pheno_startbare_ground   ***
log_sl_scaled:pheno_startdry_grassland *  
log_sl_scaled:pheno_startopen_forest   ***
log_sl_scaled:pheno_starttree_grass    ***
cos_ta_scaled:pheno_startbare_ground   ***
cos_ta_scaled:pheno_startdry_grassland    
cos_ta_scaled:pheno_startopen_forest   ***
cos_ta_scaled:pheno_starttree_grass    ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
toc()
82.413 sec elapsed
Code
beep(sound = 2)

# plot result as incidence rate ratios
# plot_model(buffalo.pheno.wet.mvmt)
plot_model(buffalo.pheno.wet.mvmt, transform = NULL)

Code
saveRDS(buffalo.pheno.wet.mvmt, file = paste0("buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_wet_", Sys.Date(), ".rds"))
paste0("Model location: buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_wet_", Sys.Date(), ".rds")
[1] "Model location: buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_wet_2025-05-20.rds"

Plot both models

Code
# plot_model(buffalo.pheno.dry.mvmt)
# plot_model(buffalo.pheno.wet.mvmt)
# 
# plot_model(buffalo.pheno.dry.mvmt, transform = NULL)
# plot_model(buffalo.pheno.wet.mvmt, transform = NULL)

Create manuscript figures

Code
coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
                      "Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass", 
                                      "sl", "log(sl)", "cos(ta)", 
                                      "sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
                                      "log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
                                      "cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
                                      ),
                      "Estimate" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"],
                      "SE" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Std. Error"]
) %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
  if (coef_df$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))

# Specify the order in which the coefficients should be plotted
order <- c(
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass",
  "sl",
  "log(sl)",
  "cos(ta)",
  "sl:Floodplain",
  "sl:DryGrassland",
  "sl:OpenForest",
  "sl:TreeGrass",
  "log(sl):Floodplain",
  "log(sl):DryGrassland",
  "log(sl):OpenForest",
  "log(sl):TreeGrass",
  "cos(ta):Floodplain",
  "cos(ta):DryGrassland",
  "cos(ta):OpenForest",
  "cos(ta):TreeGrass"
)

# Specify the order in which the coefficients should be plotted
order_no_mvmt <- c(
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass"
)

# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Plot the coefficients

Code
min_x <- -0.5
max_x <- 0.8

# coefficients on the x-axis
dry_mvmt_plot <- ggplot(data = coef_df, 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 12.5, yend   = 12.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 15.5, yend   = 15.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  # drop the filter argument when plotting all covariates
  geom_errorbarh(data = coef_df2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  # ggtitle("Dry season") +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_mvmt_plot
Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text()`).

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_dry_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 150, units = "mm", scale = 1, dpi = 600)
Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text()`).

Drop the filter argument when plotting all covariates

Code
dry_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")), 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  # drop the filter argument when plotting all covariates
  geom_errorbarh(data = coef_df2_no_mvmt,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  # scale_y_discrete(limits = rev(order)) + 
  scale_y_discrete(limits = rev(order_no_mvmt)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_mvmt_plot
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_dry_mvmt_WITHOUT_MVMT_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 60, units = "mm", scale = 1, dpi = 600)
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
                      "Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass", 
                                      "sl", "log(sl)", "cos(ta)", 
                                      "sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
                                      "log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
                                      "cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
                                      ),
                      "Estimate" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"],
                      "SE" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Std. Error"]
) %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
  if (coef_df$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))

# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Plot the coefficients

Code
# coefficients on the x-axis
wet_mvmt_plot <- ggplot(data = coef_df, 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 12.5, yend   = 12.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 15.5, yend   = 15.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  geom_point(shape = 3, size = 2.5) +
  geom_errorbarh(data = coef_df2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_mvmt_plot
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_wet_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 150, units = "mm", scale = 1, dpi = 600)
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).

Drop the filter argument when plotting all covariates

Code
wet_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")), 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  geom_errorbarh(data = coef_df2_no_mvmt,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order_no_mvmt)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_mvmt_plot
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_wet_mvmt_WITHOUT_MVMT_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 60, units = "mm", scale = 1, dpi = 600)
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Updating movement parameters, with vegetation interactions

Code
# from the step generation script
tentative_shape <- 0.438167
tentative_scale <- 534.3507
tentative_kappa <- 0.1848126

sl_extent <- 1000 # metres
sl_dry <- data.frame(x = seq(0, sl_extent, length.out = 1001))
Code
coef_df_dry_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
                           "scaled_coefficient" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"])

# step length coefficients
sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate)) 
sl_dry_mvmt_mvparams <- sl_dry_mvmt_mvparams %>% mutate(
  scale_factor = sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
  
# log step length coefficients
log_sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate)) 
log_sl_dry_mvmt_mvparams <- log_sl_dry_mvmt_mvparams %>% mutate(
  scale_factor = log_sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# turning angle coefficients
cos_ta_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate)) 
cos_ta_dry_mvmt_mvparams <- cos_ta_dry_mvmt_mvparams %>% mutate(
  scale_factor = cos_ta_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

tentative_mean_sl <- tentative_shape * tentative_scale

# shrub_grass
shrub_grass_mvmt_dry <- data.frame(
  "pheno" = "ShrubGrass",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1],
  "updated_scale" = 1/((1/tentative_scale) - sl_dry_mvmt_mvparams$nat_coef[1]),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1]
)

shrub_grass_mean_sl_dry <- shrub_grass_mvmt_dry$updated_shape * shrub_grass_mvmt_dry$updated_scale

# floodplain
floodplain_mvmt_dry <- data.frame(
  "pheno" = "Floodplain",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[2],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[2])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[2]
)

floodplain_mean_sl_dry <- floodplain_mvmt_dry$updated_shape * floodplain_mvmt_dry$updated_scale

# dry_grassland
dry_grassland_mvmt_dry <- data.frame(
  "pheno" = "DryGrassland",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[3],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[3])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[3]
)

dry_grassland_mean_sl_dry <- dry_grassland_mvmt_dry$updated_shape * dry_grassland_mvmt_dry$updated_scale

# open_forest
open_forest_mvmt_dry <- data.frame(
  "pheno" = "OpenForest",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[4],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[4])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[4]
)

open_forest_mean_sl_dry <- open_forest_mvmt_dry$updated_shape * open_forest_mvmt_dry$updated_scale

# tree_grass
tree_grass_mvmt_dry <- data.frame(
  "pheno" = "TreeGrass",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[5],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[5])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[5]
)

tree_grass_mean_sl_dry <- tree_grass_mvmt_dry$updated_shape * tree_grass_mvmt_dry$updated_scale

Creating Gamma distributions

Code
sl_dry_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))

# Tentative distribution
sl_dry_mvmt$Tentative <- dgamma(x = sl_dry_mvmt$x,
                                 shape = tentative_shape,
                                 scale = tentative_scale)
# Shrub grass
sl_dry_mvmt$ShrubGrass <- dgamma(x = sl_dry_mvmt$x,
                             shape = shrub_grass_mvmt_dry$updated_shape,
                             scale = shrub_grass_mvmt_dry$updated_scale)
# Bare ground
sl_dry_mvmt$Floodplain <- dgamma(x = sl_dry_mvmt$x,
                             shape = floodplain_mvmt_dry$updated_shape,
                             scale = floodplain_mvmt_dry$updated_scale)
# Dry grassland
sl_dry_mvmt$DryGrassland <- dgamma(x = sl_dry_mvmt$x,
                             shape = dry_grassland_mvmt_dry$updated_shape,
                             scale = dry_grassland_mvmt_dry$updated_scale)
# Open forest
sl_dry_mvmt$OpenForest <- dgamma(x = sl_dry_mvmt$x,
                             shape = open_forest_mvmt_dry$updated_shape,
                             scale = open_forest_mvmt_dry$updated_scale)
# Tree grass
sl_dry_mvmt$TreeGrass <- dgamma(x = sl_dry_mvmt$x,
                             shape = tree_grass_mvmt_dry$updated_shape,
                             scale = tree_grass_mvmt_dry$updated_scale)

# Pivot from wide to long data
sl_dry_mvmt <- sl_dry_mvmt %>% 
  pivot_longer(cols = -x)

# Specify the order in which the coefficients should be plotted
order <- c(
  "ShrubGrass",
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass",
  "Tentative"
)

sl_dry_mvmt$name <- factor(sl_dry_mvmt$name, levels = order)


mean_values <- c("ShrubGrass" = round(shrub_grass_mean_sl_dry, 1),
                 "Floodplain" = round(floodplain_mean_sl_dry, 1), 
                 "DryGrassland" = round(dry_grassland_mean_sl_dry, 1), 
                 "OpenForest" = round(open_forest_mean_sl_dry, 1), 
                 "TreeGrass" = round(tree_grass_mean_sl_dry, 1), 
                 "Tentative" = round(tentative_mean_sl, 1))

# Construct dynamic labels
labels <- paste(names(mean_values), ": Mean SL = ", mean_values, sep = "")

# Plot
ggplot(sl_dry_mvmt, aes(x = x, y = value, color = name, linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
  scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black"),
                      labels = labels) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted"),
                        labels = labels) +
  theme_bw() +
  theme(legend.position = c(0.65, 0.65))
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_sl_dry_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Creating von Mises distributions

Code
ta_dry_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))

# Tentative distribution
ta_dry_mvmt$Tentative <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = tentative_kappa)
# Shrub grass
ta_dry_mvmt$ShrubGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = shrub_grass_mvmt_dry$updated_kappa)
# Bare ground
ta_dry_mvmt$Floodplain <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = floodplain_mvmt_dry$updated_kappa)
# Dry grassland
ta_dry_mvmt$DryGrassland <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = dry_grassland_mvmt_dry$updated_kappa)
# Open forest
ta_dry_mvmt$OpenForest <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = open_forest_mvmt_dry$updated_kappa)
# Tree grass
ta_dry_mvmt$TreeGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = tree_grass_mvmt_dry$updated_kappa)

# Pivot from wide to long data
ta_dry_mvmt <- ta_dry_mvmt %>% 
  pivot_longer(cols = -x)

ta_dry_mvmt$name <- factor(ta_dry_mvmt$name, levels = order)

# Plot
ggplot(ta_dry_mvmt, 
       aes(x = x, y = value, color = factor(name), linetype = factor(name))
       ) +
  geom_line(alpha = 0.75, linewidth = 1) +
  xlab("Turning angle (radians)") +
  scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
  scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
                     labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black")) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted")) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Dry season") +
  theme_bw() +
  theme(legend.position = "bottom")

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_ta_dry_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 600)
Code
coef_df_wet_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
                                        "scaled_coefficient" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"])

# step length coefficients
sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate)) 
sl_wet_mvmt_mvparams <- sl_wet_mvmt_mvparams %>% mutate(
  scale_factor = sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# log step length coefficients
log_sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate)) 
log_sl_wet_mvmt_mvparams <- log_sl_wet_mvmt_mvparams %>% mutate(
  scale_factor = log_sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# turning angle coefficients
cos_ta_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate)) 
cos_ta_wet_mvmt_mvparams <- cos_ta_wet_mvmt_mvparams %>% mutate(
  scale_factor = cos_ta_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# shrub_grass
shrub_grass_mvmt_wet <- data.frame(
  "pheno" = "ShrubGrass",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1],
  "updated_scale" = 1/((1/tentative_scale) - sl_wet_mvmt_mvparams$nat_coef[1]),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1]
)

shrub_grass_mean_sl_wet <- shrub_grass_mvmt_wet$updated_shape * shrub_grass_mvmt_wet$updated_scale

# floodplain
floodplain_mvmt_wet <- data.frame(
  "pheno" = "Floodplain",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[2],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[2])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[2]
)

floodplain_mean_sl_wet <- floodplain_mvmt_wet$updated_shape * floodplain_mvmt_wet$updated_scale

# dry_grassland
dry_grassland_mvmt_wet <- data.frame(
  "pheno" = "wetGrassland",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[3],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[3])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[3]
)

dry_grassland_mean_sl_wet <- dry_grassland_mvmt_wet$updated_shape * dry_grassland_mvmt_wet$updated_scale

# open_forest
open_forest_mvmt_wet <- data.frame(
  "pheno" = "OpenForest",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[4],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[4])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[4]
)

open_forest_mean_sl_wet <- open_forest_mvmt_wet$updated_shape * open_forest_mvmt_wet$updated_scale

# tree_grass
tree_grass_mvmt_wet <- data.frame(
  "pheno" = "TreeGrass",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[5],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[5])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[5]
)

tree_grass_mean_sl_wet <- tree_grass_mvmt_wet$updated_shape * tree_grass_mvmt_wet$updated_scale

Creating Gamma distributions

Code
sl_wet_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))

# Tentative distribution
sl_wet_mvmt$Tentative <- dgamma(x = sl_wet_mvmt$x,
                                shape = tentative_shape,
                                scale = tentative_scale)
# Shrub grass
sl_wet_mvmt$ShrubGrass <- dgamma(x = sl_wet_mvmt$x,
                                 shape = shrub_grass_mvmt_wet$updated_shape,
                                 scale = shrub_grass_mvmt_wet$updated_scale)
# Bare ground
sl_wet_mvmt$Floodplain <- dgamma(x = sl_wet_mvmt$x,
                                 shape = floodplain_mvmt_wet$updated_shape,
                                 scale = floodplain_mvmt_wet$updated_scale)
# Dry grassland
sl_wet_mvmt$DryGrassland <- dgamma(x = sl_wet_mvmt$x,
                                   shape = dry_grassland_mvmt_wet$updated_shape,
                                   scale = dry_grassland_mvmt_wet$updated_scale)
# Open forest
sl_wet_mvmt$OpenForest <- dgamma(x = sl_wet_mvmt$x,
                                 shape = open_forest_mvmt_wet$updated_shape,
                                 scale = open_forest_mvmt_wet$updated_scale)
# Tree grass
sl_wet_mvmt$TreeGrass <- dgamma(x = sl_wet_mvmt$x,
                                shape = tree_grass_mvmt_wet$updated_shape,
                                scale = tree_grass_mvmt_wet$updated_scale)

# Pivot from wide to long data
sl_wet_mvmt <- sl_wet_mvmt %>% 
  pivot_longer(cols = -x)

sl_wet_mvmt$name <- factor(sl_wet_mvmt$name, levels = order)

mean_values_wet <- c("ShrubGrass" = round(shrub_grass_mean_sl_wet, 1),
                     "Floodplain" = round(floodplain_mean_sl_wet, 1), 
                     "DryGrassland" = round(dry_grassland_mean_sl_wet, 1), 
                     "OpenForest" = round(open_forest_mean_sl_wet, 1), 
                     "TreeGrass" = round(tree_grass_mean_sl_wet, 1), 
                     "Tentative" = round(tentative_mean_sl, 1))

# Construct dynamic labels
labels <- paste(names(mean_values_wet), ": Mean SL = ", mean_values_wet, sep = "")

# Plot
ggplot(sl_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
  scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black"),
                      labels = labels) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted"),
                        labels = labels) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Wet season") +
  theme_bw() +
  theme(legend.position = c(0.65, 0.65))

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_sl_wet_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Creating von Mises distributions

Code
ta_wet_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))

# Tentative distribution
ta_wet_mvmt$Tentative <- circular::dvonmises(x = ta_wet_mvmt$x,
                                             mu = 0,
                                             kappa = tentative_kappa)
# Shrub grass
ta_wet_mvmt$ShrubGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = shrub_grass_mvmt_wet$updated_kappa)
# Bare ground
ta_wet_mvmt$Floodplain <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = floodplain_mvmt_wet$updated_kappa)
# Dry grassland
ta_wet_mvmt$DryGrassland <- circular::dvonmises(x = ta_wet_mvmt$x,
                                                mu = 0,
                                                kappa = dry_grassland_mvmt_wet$updated_kappa)
# Open forest
ta_wet_mvmt$OpenForest <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = open_forest_mvmt_wet$updated_kappa)
# Tree grass
ta_wet_mvmt$TreeGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
                                             mu = 0,
                                             kappa = tree_grass_mvmt_wet$updated_kappa)

# Pivot from wide to long data
ta_wet_mvmt <- ta_wet_mvmt %>% 
  pivot_longer(cols = -x)

ta_wet_mvmt$name <- factor(ta_wet_mvmt$name, levels = order)

# Plot
ggplot(ta_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  xlab("Turning angle (radians)") +
  scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
  scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
                     labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black")) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted")) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Dry season") +
  theme_bw() +
  theme(legend.position = "bottom")

Code
ggsave(paste0("outputs/plots/pheno_model_mixed_diag_ta_wet_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 600)

Session Information

Code
sessionInfo()
R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Brisbane
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] amt_0.2.2.0         circular_0.5-1      lemon_0.5.0        
 [4] beepr_2.0           tictoc_1.2.1        glmmTMB_1.1.11     
 [7] sjPlot_2.8.17       TwoStepCLogit_1.2.6 lubridate_1.9.4    
[10] forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4        
[13] purrr_1.0.4         readr_2.1.5         tidyr_1.3.1        
[16] tibble_3.2.1        ggplot2_3.5.2       tidyverse_2.0.0    

loaded via a namespace (and not attached):
 [1] sjlabelled_1.2.0    tidyselect_1.2.1    farver_2.1.2       
 [4] fastmap_1.2.0       bayestestR_0.15.3   sjstats_0.19.0     
 [7] digest_0.6.37       timechange_0.3.0    lifecycle_1.0.4    
[10] sf_1.0-20           survival_3.8-3      magrittr_2.0.3     
[13] compiler_4.5.0      rlang_1.1.6         tools_4.5.0        
[16] yaml_2.3.10         knitr_1.50          labeling_0.4.3     
[19] bit_4.6.0           classInt_0.4-11     plyr_1.8.9         
[22] RColorBrewer_1.1-3  KernSmooth_2.23-26  withr_3.0.2        
[25] numDeriv_2016.8-1.1 grid_4.5.0          datawizard_1.1.0   
[28] e1071_1.7-16        scales_1.4.0        MASS_7.3-65        
[31] insight_1.2.0       cli_3.6.5           mvtnorm_1.3-3      
[34] crayon_1.5.3        rmarkdown_2.29      ragg_1.4.0         
[37] reformulas_0.4.1    generics_0.1.3      rstudioapi_0.17.1  
[40] performance_0.13.0  tzdb_0.5.0          parameters_0.25.0  
[43] minqa_1.2.8         DBI_1.2.3           proxy_0.4-27       
[46] audio_0.1-11        splines_4.5.0       parallel_4.5.0     
[49] effectsize_1.0.0    vctrs_0.6.5         boot_1.3-31        
[52] Matrix_1.7-3        jsonlite_2.0.0      hms_1.1.3          
[55] bit64_4.6.0-1       systemfonts_1.2.3   units_0.8-7        
[58] glue_1.8.0          nloptr_2.2.1        stringi_1.8.7      
[61] gtable_0.3.6        ggeffects_2.2.1     lme4_1.1-37        
[64] pillar_1.10.2       htmltools_0.5.8.1   R6_2.6.1           
[67] TMB_1.9.17          textshaping_1.0.1   Rdpack_2.6.4       
[70] vroom_1.6.5         evaluate_1.0.3      lattice_0.22-6     
[73] rbibutils_2.3       class_7.3-23        Rcpp_1.0.14        
[76] gridExtra_2.3       nlme_3.1-168        mgcv_1.9-1         
[79] xfun_0.52           sjmisc_2.8.10       pkgconfig_2.0.3